import string
import scipy
import Tkinter, tkFileDialog
import numpy as np
import pandas as pd
from scipy.spatial.distance import cdist
from scipy import ndimage
import skimage
import skimage.filters
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import matplotlib.image as mpimg
from PIL import Image
import cmath
import os
import pdb
import sys
import glob
sys.path.append(os.path.abspath("C:\Users\Scherer Lab E\Documents\GitHub\Python_Data_Analysis"))
from common_functions import *
%matplotlib inline

Define Useful Functions

def makeGaussian(size, fwhm = 3, center=None):
    """ Make a square gaussian kernel.
    size is the length of a side of the square
    fwhm is full-width-half-maximum, which
    can be thought of as an effective radius.
    """

    x = np.arange(0, size, 1, float)
    y = x[:,np.newaxis]
    
    if center is None:
        x0 = y0 = size // 2
    else:
        x0 = center[0]
        y0 = center[1]
    
    return np.exp(-4*np.log(2) * ((x-x0)**2 + (y-y0)**2) / fwhm**2)

def normalize_array(array):
    return array/np.sum(array)

def make_complex_phase_amp(amp, phase):
    return amp * (np.cos(phase) + 1j*np.sin(phase))

Algorithm

def aa_array_generation(beamlet_array_pts, spacing, phase_mask_size, additive_term=0.5, wavelength=800):
    '''
    This function uses an Additive Adaptive algorithm to compute the phase mask
    to do an array of beamlets about the center of the mask. The source for this
    method is:
    
    http://physics.nyu.edu/grierlab/cgh3c/
    
    However, this algorithm is not optimized because it should only take ~8
    iterations to reach convergence but still produces a viable diffraction
    grating.
    
    :param (tuple) beamlet_array_pts: Describes the beamlet array you want to make of beamlets. The beamlets
    will automatically be around the center
    :param spacing: the spacing between each beamlet (in pixels)
    :param additive_term: The amount of our target phase that gets added to the non-target phase. Higher means
    that the target phase is used more
    :param wavelength: The wavelength of your laser light
    '''
    k = 2*np.pi /(wavelength)

    # Create the Target Array
    target_array = np.zeros([phase_mask_size,phase_mask_size])
    x_coors, y_coors = np.meshgrid((np.arange(0,beamlet_array_pts[0])*spacing-int(0.5*(beamlet_array_pts[0]-1)*spacing)),
                               (np.arange(0,beamlet_array_pts[1])*spacing-int(0.5*(beamlet_array_pts[1]-1)*spacing)))
    x_coors += phase_mask_size/2
    y_coors += phase_mask_size/2
    #pdb.set_trace()
    target_array[x_coors.flatten(), y_coors.flatten()] = 1
    
    # Use a Guassian Kernel (fwh 3/4s of phase mask size) and random phase
    # to produce the initial electric field
    init_phase = (-np.pi - np.pi) * np.random.rand(phase_mask_size, phase_mask_size) + np.pi
    init_amp = makeGaussian(phase_mask_size,phase_mask_size*3/4.0)
    init_complex = init_amp * (np.cos(init_phase) + 1j*np.sin(init_phase))
    image_plane = np.fft.fft2(make_complex_phase_amp(init_amp, init_phase))
    
    # Set initial chi and error before loop
    error_prev = None
    chi = 10**10
    
    while chi > 10**-6:
        
        # Extract Amp and Phase from Image Plane
        image_plane_amp = np.fft.fftshift(abs(image_plane))
        image_plane_phase = np.angle(image_plane)
        
        # Calculate Error and check if complete
        error_curr = (np.sum(target_array - 
                             image_plane_amp)**2)/phase_mask_size**2
        if error_prev == None:
            error_prev = error_curr
        else:
            chi = abs(error_curr - error_prev)/error_curr
            #print chi
            error_prev = error_curr
            if chi < 10**(-6):
                break
        
        # Mix the image plane amplitude with the target amplitude
        new_image_plane_amp = additive_term*normalize_array(target_array) + (1-additive_term)*normalize_array(image_plane_amp)
        
        # Compute the E Field in the SLM Plane
        slm_plane = np.fft.ifft2(make_complex_phase_amp(np.fft.fftshift(new_image_plane_amp), image_plane_phase))
        slm_plane_amp = np.abs(slm_plane)
        slm_plane_phase = np.angle(slm_plane)
        
        # Compute the E Field in the Image plane, replace amplitude in SLM plane
        # with gaussian
        image_plane = np.fft.fft2(make_complex_phase_amp(init_amp, slm_plane_phase))
    
    plt.figure(figsize=(10,10))
    plt.subplot(221)
    plt.title('Image Plane Amp')
    plt.imshow(np.abs(image_plane_amp), cmap='gray', interpolation='none')
    plt.gca().invert_yaxis()
    plt.subplot(222)
    plt.title('SLM Plane Amp')
    plt.imshow(slm_plane_amp, cmap='gray', interpolation='none')
    plt.subplot(223)
    plt.title('Image Plane Phase')
    plt.imshow(image_plane_phase, cmap='gray', interpolation='none')
    plt.subplot(224)
    plt.title('SLM Plane Phase')
    plt.imshow(slm_plane_phase, cmap='gray', interpolation='none')
    plt.tight_layout()
    plt.show()
    return slm_plane_phase, image_plane_amp
    

phase_output, img_plane_amp = aa_array_generation([1,2], 15, 600)

phase_output = phase_output+abs(np.pi)
phase_output *= 255/(2*np.pi)
phase_output = np.tile(phase_output,(2,2))[:600,:792]

plt.imshow(np.tile(phase_output,(2,2))[:600,:792])

<matplotlib.image.AxesImage at 0x1b249d68>

def wrap2value(array, lower_value, upper_value):
    """
    Wraps the phase around 2pi such that multiples of 2pi greater than 2pi will equal 2pi and multiples of
     2p less than 0 will equal 0. This is called automatically on the phase attribute
    """
    array[array % (upper_value) != lower_value] %= (upper_value)
    array[np.logical_and(array % (upper_value) == lower_value, array / (upper_value) > lower_value)] = upper_value
    array[np.logical_and(array % (upper_value) == lower_value, array / (upper_value) < lower_value)] = lower_value
    return array

save_dir = "K:\Pat's_Projects\Custom Phase Mask\Phase Masks for Delphine\\"
phase_mask = Image.fromarray(phase_output.astype('uint8'))
phase_mask.save(save_dir+'1x2sp15_diffraction_array.bmp', format='bmp')